#pragma once

void ludcmp(complex **a, int n, int *indx,int iflag)
{
	double tiny=1e-20;

	complex *vv=(complex*)malloc(sizeof(complex)*(n));
	complex sum;
	double aamax;
	complex dum;
        int imax=0;	
	for (int i=0; i<n; i++)
	{
		aamax=0.;
		for (int j=0; j<n; j++)
		{
			if (cabs(a[i][j])>aamax) aamax=cabs(a[i][j]);

		}

		if (cabs(aamax)<0.0001) { iflag=1; break;}
		vv[i]=1./aamax; 
	}
	for (int j=0; j<n; j++)
	{
		for(int i=0; i<j; i++)
		{
			sum=a[i][j];
			for (int k=0; k<i; k++) sum=sum-a[i][k]*a[k][j];
			a[i][j]=sum;
		}
		aamax=0.;	
	      	for (int i=j; i<n; i++)
		{
			 sum=a[i][j];
                         for (int k=0; k<j; k++) sum=sum-a[i][k]*a[k][j];
                         a[i][j]=sum;

			dum=vv[i]*cabs(sum);
			if (cabs(dum)>=aamax)
			{
					imax=i;
					aamax=dum;
			}

		}
		if (j != imax+1)
		{

			for (int k=0; k<n; k++)
			{
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (j !=(n-1)) 
		{
			if (a[j][j]==0.0 ) a[j][j]=tiny;
			dum=1./a[j][j];
			for (int i=j+1; i<n; i++) a[i][j]=a[i][j]*dum;			

		}

	}
	if (a[n-1][n-1]==0.0)  a[n-1][n-1]=tiny;

	
}


void lubksb(complex **a, int n,int *indx, complex *b)

{
	complex sum;

	int ll;
	
	int ii=-1;
	for (int i=0; i<n; i++)
	{
		ll=indx[i];
		sum=b[ll]; 
		b[ll]=b[i];
		if (ii != 1)
		{
			for (int j=ii; j<i; j++) sum=sum-a[i][j]*b[j];
		}
		else
		{
			if (sum != 0.0) ii=i;
		}
		b[i]=sum;
	}

	for (int i=(n-1); i>=0 ; i--)
	{
		sum=b[i]; 
		for (int j=i+1; j<n; j++) sum=sum-a[i][j]*b[j];
		b[i]=sum/a[i][i];

	}
}

void inverse(complex **a, int n,complex **inva,int iflag)
{
//     From Numerical Recipes (page 38)

//     Input Parameters:   a  = matrix to be inverted
//                         n  = order of the matrix a
//
//     Output Parameters:  inva = inverse of the matrix a
//                         a = original matrix
//                             (ie. it's not overwritten)
//                         iflag = 0 if inverse found, =1 if a is singular



      int indx[n];

      complex **aa=(complex**)malloc(sizeof(complex)*(n));
      for(int i=0; i<n; i++) aa[i]=(complex*)malloc(sizeof(complex)*(n));

      complex **iinva=(complex**)malloc(sizeof(complex)*(n));
      for(int i=0; i<n; i++) iinva[i]=(complex*)malloc(sizeof(complex)*(n));

      complex *col=(complex*)malloc(sizeof(complex)*(n));


      iflag=0; //reset flag

       //Reescale
	double max;
	max=cabs(a[0][0]);

	for (int i=0; i<n; i++)
	{
		for (int j=0; j<n; j++)
		{
			iinva[i][j]=0.;
			aa[i][j]=a[i][j];
			if (cabs(a[i][j]) > max) max=cabs(a[i][j]);
		}
		iinva[i][i]=1.;
	}
	for (int i=0; i<n; i++) for (int j=0; j<n; j++) aa[i][j]=aa[i][j]/max;

      //Decompose the matrix
      ludcmp(aa,n,indx,iflag);
      for(int j=0;j<n;j++)
    	{
	  //find inverse by columns
    	  for(int i=0;i<n;i++) col[i]=0.;
     	  col[j]=1.0;
     	  lubksb(aa,n,indx,col);
	  for(int i=0;i<n;i++) iinva[i][j]=col[i];
        }
       for (int i=0; i<n; i++) for (int j=0; j<n; j++) inva[i][j]=iinva[i][j]/max;


}
////////////////////////////////////////////////
//REAL	
//////////////////////////////////////////////


void ludcmp_real(double **a, int n, int *indx,int iflag)
{
        double tiny=1e-20;

        double *vv=(double*)malloc(sizeof(double)*(n));
        double sum;
        double aamax;
        double dum;
        int imax=0;
        for (int i=0; i<n; i++)
        {
                aamax=0.;
                for (int j=0; j<n; j++)
                {
                        if (fabs(a[i][j])>aamax) aamax=fabs(a[i][j]);

                }

                if (fabs(aamax)<0.0001) { iflag=1; break;}
                vv[i]=1./aamax;
        }
        for (int j=0; j<n; j++)
        {
                for(int i=0; i<j; i++)
                {
                        sum=a[i][j];
                        for (int k=0; k<i; k++) sum=sum-a[i][k]*a[k][j];
                        a[i][j]=sum;
                }
                aamax=0.;
                for (int i=j; i<n; i++)
                {
                         sum=a[i][j];
                         for (int k=0; k<j; k++) sum=sum-a[i][k]*a[k][j];
                         a[i][j]=sum;

                        dum=vv[i]*fabs(sum);
                        if (fabs(dum)>=aamax)
                        {
                                        imax=i;
                                        aamax=dum;
                        }

                }
                if (j != imax+1)
                {

                        for (int k=0; k<n; k++)
                        {
                                dum=a[imax][k];
                                a[imax][k]=a[j][k];
                                a[j][k]=dum;
                        }
                        vv[imax]=vv[j];
                }
                indx[j]=imax;
                if (j !=(n-1))
                {
                        if (a[j][j]==0.0 ) a[j][j]=tiny;
                        dum=1./a[j][j];
                        for (int i=j+1; i<n; i++) a[i][j]=a[i][j]*dum;

                }

        }
        if (a[n-1][n-1]==0.0)  a[n-1][n-1]=tiny;


}


void lubksb_real(double **a, int n,int *indx, double *b)

{
        double sum;

        int ll;

        int ii=-1;
        for (int i=0; i<n; i++)
        {
                ll=indx[i];
                sum=b[ll];
                b[ll]=b[i];
                if (ii != 1)
                {
                        for (int j=ii; j<i; j++) sum=sum-a[i][j]*b[j];
                }
                else
                {
                        if (sum != 0.0) ii=i;
                }
                b[i]=sum;
        }

        for (int i=(n-1); i>=0 ; i--)
        {
                sum=b[i];
                for (int j=i+1; j<n; j++) sum=sum-a[i][j]*b[j];
                b[i]=sum/a[i][i];

        }
}

void inverse_real(double **a, int n,double **inva,int iflag)
{
//     From Numerical Recipes (page 38)

//     Input Parameters:   a  = matrix to be inverted
//                         n  = order of the matrix a
//
//     Output Parameters:  inva = inverse of the matrix a
//                         a = original matrix
//                             (ie. it's not overwritten)
//                         iflag = 0 if inverse found, =1 if a is singular



      int indx[n];

      double **aa=(double**)malloc(sizeof(double)*(n));
      for(int i=0; i<n; i++) aa[i]=(double*)malloc(sizeof(double)*(n));

      double **iinva=(double**)malloc(sizeof(double)*(n));
      for(int i=0; i<n; i++) iinva[i]=(double*)malloc(sizeof(double)*(n));

      double *col=(double*)malloc(sizeof(double)*(n));


      iflag=0; //reset flag

       //Reescale
        double max;
        max=fabs(a[0][0]);

        for (int i=0; i<n; i++)
        {
                for (int j=0; j<n; j++)
                {
                        iinva[i][j]=0.;
                        aa[i][j]=a[i][j];
                        if (fabs(a[i][j]) > max) max=fabs(a[i][j]);
                }
                iinva[i][i]=1.;
        }
        for (int i=0; i<n; i++) for (int j=0; j<n; j++) aa[i][j]=aa[i][j]/max;

      //Decompose the matrix
      ludcmp_real(aa,n,indx,iflag);
      for(int j=0;j<n;j++)
        {
          //find inverse by columns
          for(int i=0;i<n;i++) col[i]=0.;
          col[j]=1.0;
          lubksb_real(aa,n,indx,col);
          for(int i=0;i<n;i++) iinva[i][j]=col[i];
        }
       for (int i=0; i<n; i++) for (int j=0; j<n; j++) inva[i][j]=iinva[i][j]/max;


}

